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A model for ionic solutions with an attractive short-range pair interaction between the ions is 
presented. The short-range interaction is accounted for by adding a quadratic non-local term to the 
Poisson-Boltzmann free energy. The model is used to study solvent effects in a planar electric double 
layer. The counter-ion density is found to increase near the charged surface, as compared with the 
Poisson-Boltzmann theory, and to decrease at larger distances. The ion density profile is studied 
analytically in the case where the ion distribution near the plate is dominated only by counter-ions. 
Further away from the plate the density distribution can be described using a Poisson-Boltzmann 
theory with an effective surface charge that is smaller than the actual one. 
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I. INTRODUCTION 



Electrolytes, in contact with charged surfaces or 
macro-ions, play an important role in determining the 
properties of many biological and chemical systems. One 
of the most widely used tools for studying ions in aque- 
ous solutions is the Poisson-Boltzmann (PB) theory [Q-l). 
The mathematical and conceptual simplicity of this the- 
ory makes it very appealing both for numerical computa- 
tion [Q , and for gaining insight into the underlying phys- 
ical principles. Although the theory contains important 
simplifications, it has proven to be a useful and accurate 
tool in the study of systems such as colloidal dispersions 
[H,D, biological membranes [D, polyelectrolytes |Q and 
complex systems, e.g., proteins or DNA interacting with 
charged membranes |ll| |l3| . 

The Poisson-Boltzmann theory is obtained by making 
two simplifying approximations. The first approxima- 
tion is the treatment of the electrostatic interactions on 
a mean-field level. The ions are treated as independent 
charged particles interacting with an external electro- 
static potential, derived self-consistently from the mean 
charge density distribution. Thus, correlations between 
the ion positions are not taken into account. The second 
approximation is the treatment of the ions as point-like 
objects, interacting only through the electrostatic inter- 
action in a dielectric medium. In reality, ions in aque- 
ous solutions have more intricate interactions These 
include a non-Coulombic interaction between ion pairs, 
which is mainly a short-range steric repulsion, interac- 
tions with the polar solvent molecules and short-range 
interactions with the confining charged surfaces. 



Various models have been proposed for the inclusion 
of effects not accounted for by the PB theory. These 
include liquid state theory approaches |r^-[T^, field the- 
ory expansions (l^,^^ , computer simulatio nspi] - ^ and 
other modifications to the PB theory p^27[|. Most of 
these models remain within the framework of the so- 
called "primitive model" , in which the interaction be- 
tween the ions is modeled as a purely repulsive hard- 
core interaction. On the other hand, relatively few works 
have addressed explicitly the discrete nature of the sol- 
vent molecules p^ , p8| -p0[. Clearly, the replacement of 
the solvent by a continuous medium cannot be precise 
when the inter-ion distance is comparable to the solvent 
molecular size. Therefore, when the ions reach high den- 
sities the discreteness of the solvent is expected to have 
an important effect on the ionic distribution. This is of 
particular importance for water. Due to its high polarity, 
the strong screening of the electrostatic interaction (rep- 
resented by the dielectric constant) is modified at small 
ion separations. 

Using the surface force apparatus ||], it is possible 
to measure precisely the force between charged mica 
plates. These measurements supply evidence for the im- 
portance of the solvent structure in aqueous solutions 
|3l| , p2| . At inter-plate separations below approximately 
20 A^ignificant deviations are found from the prediction 
of the Derjaguin-Landau-Verwey-Overbeek (DLVO) the- 
ory p|j3^ . The measured force is oscillatory or consists of 
a series of steps, with a period corresponding to the wa- 
ter molecular size. Oscillatory forces are known to arise 
as a result of the solvent structuring in layers between 
surfaces However, a repulsive contribution is found 
in addition to the oscillatory force at plate separations 
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below several nanometers pi],B2| . This repulsive force is 



often referred to as the "hydration force" 
origin is not yet completely understood \U9 




and its 



A. Aqueous pair potential model 

Recently |^,^, an aqueous pair potential model has 
been proposed for electrolytes, in which the effect of the 
solvent on the ions is described as a short-range two-body 
interaction between the ions. The solvent is replaced by 
a continuum dielectric medium as in PB theory, but the 
ions also interact through a two-body short-range hy- 
dration interaction Q]. This is shown schematically in 
Fig. 1. 

This aqueous pair potential model , involves sev- 

eral simplifying assumptions. One is that the effect of the 
solvent can be represented as a linear superposition of 
two-body potentials between all ion pairs. Another sim- 
plification is that the effective potential between the ions 
is taken as the effective potential in the bulk, regardless 
of the ion concentration, and of the geometry imposed by 
the charged surfaces. Finally, a short-range surface-ion 
effective potential should be included in addition to the 
ion-ion effective potential. Despite of the simplifications 
made in the aqueous pair potential model, it offers a first 
step towards a qualitative understanding of solvent ef- 
fects on the ion distribution, in particular near highly 
charged surfaces. 



B. Effective ion pair interaction 

For the short-range ion-ion interaction, the so-called 
potential of mean force between ions in solution can 
be used. Potentials of mean force are defined as 
—ksT log gij{r) where (7y(r) are the ion-ion radial dis- 
tribution functions for ion pairs of species i and j . The 
radial distribution functions have been calculated numer- 
ically for a single ion pair immersed in an aqueous solu- 
tion using molecular dynamics techniques [p6| - |38[ . 

An alternative approach has been proposed in Refs. 
p9| , ^ . In this approach, a Hamiltonian consisting of a 
pairwise effective potential between the ions is obtained 
using the so-called "reverse Monte-Carlo" approach. The 
ion-ion radial distribution functions are first calculated 
using a molecular dynamics simulation for a system in- 
cluding solvent molecules and a finite concentration of 
ions. The ion-ion effective potential in the system with- 
out the solvent is then adjusted iteratively until the same 
distribution functions are obtained using Monte-Carlo 
simulations. 

The different available calculations of potentials of 
mean force differ in their quantitative predictions. This 
may be a result of high sensitivity of the models to de- 
tailed features used for the water molecules and for the 
inter-molecular interactions [^ . However, all the po- 
tentials of mean force as well as the effective potentials 
p9| are qualitatively similar j4|]. Thus, for the purpose 
of the present work, aiming at a qualitative understand- 
ing of solvent effects, any one of these potentials may be 
used. 




FIG. 1. Schematic description of the aqueous pair potential model. An aqueous ionic solution in contact with a charged 
plate in (a) is replaced in (b) by ions in a continuum dielectric medium having a dielectric constant e, with electrostatic 
and short -range interactions Uij(r) — iiij(|r|). The z coordinate designates the distance from the charged plate, with z — 
corresponding to the distance of closest approach of the ions to the plate. The distance of closest approach is equal to dhc/2, 
where dhc is the hard-core diameter of the ions. 
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FIG. 2. Short-range effective potential between Na ion 
pairs, adapted from Ref. using simulations in a bulk NaCl 
solution of concentration 0.55 M, at room temperature |Q. 
The potential is shown in units of ksT, as a function of the 
distance between the ion centers. For ion separations smaller 
than 2.9 A a hard core interaction was assumed. The Coulomb 
interaction is subtracted to show only the short-range hydra- 
tion effect due to the water molecules. 

At large ionic separations the ion-ion effective poten- 
tial is well approximated by a screened electrostatic in- 
teraction, with the water dielectric constant in the con- 
tinuum limit. At short ionic separation, the difference 
between the total effective potential and the screened 
electrostatic interaction is a short-range potential reflect- 
ing the structure of water molecules in the ion vicinity. 
Fig. 2 shows the short-range contribution (excluding the 
screened electrostatic part) to the efi^ective potential cal- 
culated between Na^ - Na^ pairs in the reverse Monte- 
Carlo approach Below about 3 A, the electrostatic 
repulsion between the ions becomes unscreened. There- 
fore it is much larger than the screened repulsion in the 
dielectric medium, and the effective potential is strongly 
repulsive. The unscreened electrostatic potential leads to 
an effectively enlarged hard-core separation between the 
ions, relative to a hard-core diameter of about 2.3 A used 
in the short-range part of the bare ion-ion potential. At 
larger separations, the effective potential is oscillatory, 
and mainly attractive. It has a distinct minimum at an 
ion-ion separation of about 3.6 A, followed by a maximum 
and a second minimum at approximately 6 A. 



C. The present work 

The replacement of the discrete solvent by a continuum 
medium, with electrostatic and short range interactions 
between the ions, is a considerable simplification. Still, 
the statistical mechanical treatment of an electrolyte so- 
lution in this model is difficult, and requires the use of 
further approximations, or simulations. 

The Anisotropic Hyper-Netted Chain approximation 



(AHNC) Q was previously used to calculate the ef- 
fects of hydration interactions in the aqueous pair poten- 
tial model |4^,^^. When the ion concentration is large 
enough, e.g., near a highly charged surface, the hydra- 
tion interaction is found to have a significant effect on 
the distribution of ions in the solution. It was also pro- 
posed that the so-called repulsive "hydration forces" be- 
tween surfaces arise from the ionic structure near highly 
charged surfaces. According to this description, at large 
distances from the plate, the ion distribution follows a 
PB profile with a reduced effective surface charge. When 
two plates approach each other, the ions near the sur- 
faces come into contact giving rise to an apparent new 
repulsive force. 

In the present work a simple description for ions inter- 
acting through electrostatic and short-range attractive 
interactions as mediated by the solvent molecules is in- 
troduced. We apply this description to the aqueous pair 
potential model. Our aim is limited to describe the im- 
portant effects of the short-range interaction, and not to 
provide an accurate tool for their calculation. Therefore 
our model follows the PB theory as closely as possible, 
and describes the short-range interaction using a sim- 
plified term added to the free energy. The advantage 
of this approach over more elaborate treatments such as 
the AHNC @,|2|, is that it provides relatively simple 
equations that can be treated numerically and analyti- 
cally with relative ease as well as allowing extensions to 
non-planar geometries. In a planar geometry, we show 
that the effect of the ion-ion hydration interaction can 
be understood as a perturbation over the PB results. 
An increase in the concentration of counter-ions near the 
charged surface is found, and it results in an apparent 
surface charge which is reduced relatively to the PB the- 
ory. 

In addition to the calculation of charge distributions, 
the effect of the hydration interaction on the force be- 
tween charged particles or surfaces can be studied and 
will be presented elsewhere [Q. Since we will not dis- 
cuss inter-surface forces in this paper, it is worthwhile 
to mention that the results obtained using our model 
are in good agreement with AHNC calculations |Q and 
may provide an explanation to the hydration forces as 
observed in surface force measurements ||3l| , |3^ . For high 
surface charge and plate separations up to approximately 
20 A, important modifications of the PB predictions are 
found. 

The outline of the paper is as follows. Section || 
presents the model. In section HI we apply the model 
to a single charged plate, present numerical results for 
the ion density profile and discuss the modifications to 
the PB theory due to the addition of short-range inter- 
actions. In section tV we present analytical results in 
the low salt limit. We calculate the effective PB surface 
charge and the effect of the hydration interaction on the 
density profile of counter-ions in a system with no added 
salt. 
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II. THE MODEL 

A. Free energy 

We start from an approximated free energy, written as 
a functional of the various ion densities. We choose the 
electrostatic boundary conditions to be of fixed surface 
charge densities and write the free energy as a sum of the 
usual PB term and a correction term, due to hydration, 
as will be explained below: 



(1) 



We discuss first how the PB free energy is obtained, and 
then generalize this result to include the short-range hy- 
dration interaction. 



1. Poisson-Boltzmann free energy 
The Hamiltonian of the system is: 
^ = ^ V / d^r e,p,(r)<^(r) + \ i d^r, CT(r,)0(r,) 

(2) 

where V is the the volume occupied by the electrolyte 
solution, cr(rs) is the surface charge density of immobile 
charges on the boundaries dV ^ and is the charge of the 
ith ion species. The ion densities Pi(r) are: 



(3) 



where r:^ is the position of the jth ion of the ith species, 
and the electrostatic potential (/)(r) is a function the dif- 
ferent ion positions: 



dv £ |r - i-s 



(4) 



where e = 78 is the dielectric constant of water. The 
PB theory is obtained by using a mean-field approxima- 
tion for the electrostatic interaction. The Hamiltonian 
is first replaced by a mean-field Hamiltonian, where 
the electrostatic potential (^{v) is replaced by an external 
field ^'(r). In the thermodynamic limit the free energy 
can then be written as a functional of the mean densities 
of the ion species, Ci(r) = (pi(r))MF, as follows: 



J7mf — ^TsT 



log^-1 



E 

i 



d-^r 



av 



(5) 



where k^T is the thermal energy, and C,i is the fugacity 
of the ith ion species. The mean-field approximation is 
obtained by requiring that the external potential vE'(r) is 
the thermodynamical average of (Q) in the system with 
the mean-field Hamiltonian, ^'(r) = ((/i(r))MF, *-e-: 



*(r) = ^ / 



e r — r 



d^r.^(l^ (6) 
dv e|r-rs| 



This relation is equivalent to the Poisson equation: 

Ait 



V^f = > e,c, 

e 



(7) 



supplemented by the boundary condition: 



47r 

V^* ■ n|j,^ = ""(rs) on the charged surfaces (8) 

where the normal vector fi points away from the charged 
surfaces into the volume occupied by the ionic solution. 
Using this boundary condition and Eq. (|^), the second 
and third terms of Eq. (|^) can be re-expressed as: 

= f / iy^fA^ (9) 

87r Jv 

Substituting this relation in Eq. (^) we obtain the PB 
free energy: 



^PB - ^ / (Vvl/)2d3r 



log - 1 



d-^r 



A(r) V^* 



47r 



E 



d^r 



(10) 



The first term in flpB is the electrostatic free energy and 
the second term is the entropy of the ions. The fugac- 
ity Ci, in the second term, is equal in PB theory to the 
bulk concentration Cb,i of the ith ion species, Q = c^^i, 
as for an ideal gas. For more generalized free energies, a 
different relation may exist between the fugacity of each 
ion species and its respective bulk concentration. The 
electrostatic potential is a functional of the ion den- 
sities Ci, and is determined by the Poisson equation (Q) 
and the boundary conditions (|) imposed by the surface 
charges. Alternatively, in Eq. (10|) ^' is regarded as an in- 
dependent field and a third term containing a Lagrange 
multiplier A(r) is added to Ops- The PB equilibrium 
mean densities Ci(r) result from minimizing fipB. With 
the introduction of A(r) the minimization is equivalent to 
requiring an extremum of rips with respect to the three 
fields Ci, ^' and A, subject to the boundary condition (||). 
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By requiring first an extremum of flpB with respect to ^' 
and Ci the foUowing relations are obtained: 



A = — * 

An 



and: 



a = Ci exp {-(3e,^) 



(11) 



(12) 



Where /3 — l/ik^T) and Ci = Cb,i. The extremum con- 
dition with respect to A gives the Poisson equation: 



e^Ci 



(13) 



Combining these relations we obtain the PB equation: 

e.e-'^^'* (14) 



Alternatively, the first two relations, obtained from the 
extremum condition with respect to 5" and q, can be 
substituted into Eq. (p^). Formally, this gives fipB as a 
functional of A. Using Eq. ([Ill), the expression obtained 
for rjpB can be written as a functional of 5*: 



i7pR — 



Jv JdV 

-k^T I d^r^C^e-'^^"* (15) 

where the second integration is over the charged surfaces. 
Requiring an extremum of this functional with respect to 
^! is another way to obtain the Poisson-Boltzmann equa- 
tion (0). 

A more formal derivation of the mean field, PB free 
energy, and a discussion on its generalization to sys- 
tems with non-electrostatic interactions is presented in 
Ref. Q . The PB free energy (|l^) can also be derived by 
formulating the problem using field theory methods. In 
this approach the mean-field approximation is obtained 
as the saddle point of the functional integral, and cor- 
rections due to ion-ion correlations can be obtained in a 
systematic expansion |l^,p9| . 



2. Inclusion of the hydration interaction 

As discussed in the introduction, our starting point is 
a model in which the hydration interaction, arising from 
solvent effects, is described as an effective ion-pair in- 
teraction. We denote this short-range potential between 
ions of species i and j at distance r (r). The poten- 
tial is taken as the short-range effective potential between 
ions immersed in a bulk ionic solution having a specific, 
constant concentration. Therefore, Uij(r) is assumed to 
be isotropic and does not depend on the ion positions or 
the confining geometry. 



Our aim is to treat the long range electrostatic inter- 
action on the mean-field level, as in PB theory. Thus, we 
begin by considering the free energy of a system placed 
in some arbitrary field ^'(r), where the ions interact with 
each other only through the two-body potential (r) . 
Due to the short-range nature of the hydration interac- 
tion, the free energy can be obtained from a virial expan- 
sion of the grand canonical partition function. Since we 
will be interested in highly inhomogeneous systems, we 
perform an expansion in the inhomogeneous ion density. 
The derivation is given in Appendix |^. Including terms 
up to the quadratic order in the expansion we obtain: 



log- 



ksT y y^c, 

i 

+ J ^e.Ci^d^r 

J2 1 c.(r)[/,, (r-r')c, (r')d3rdV (16) 



keT 



where *I'(r) is an external field, coupled to the ith ion 
charge density e^Ci. The short-range weighted potential 



Uij in the third term of Qh is defined as: 



-/3«i,(|r-r'|) 



(17) 



where Uij is the nominal short-range interaction potential 
between ions of species i and j. This form of describing 
the short-range interaction is a rather crude approxima- 
tion, valid only in the low density limit. Its advantage 
is its simplicity. The free energy Hh amounts to setting 
the direct correlation function C2(|r — r'|) to be equal 
to —U{\r — r'l), and all higher order direct correlation 
functions to zero p6| |. 

Having found the hydration free energy llh, the elec- 
trostatic interaction can be treated on the mean-field 
level. This is done by considering ^'(r) as the electro- 
static potential and imposing the self-consistency require- 
ment of the Poisson equation (|^). This is essentially the 
approximation we used to derive the PB equation ([l4|), 
with the difference that the free energy of a dilute, non- 
interacting ion distribution is replaced by the free energy 
fih of Eq. (p^. The result is the free energy of Eq. (|^), 
with Ail defined as follows: 



An = 



knT- 



;y'c.(r)f/,,(r-r')c,(r')d3rdV 



(18) 



We conclude this section with some remarks on the 
approach presented above. Important solvent effects are 
already introduced in the PB theory by using an elec- 
trostatic interaction with a dielectric constant £ = 78 of 
water, instead of the bare electrostatic interaction. In 
the modified model a more precise effective potential be- 
tween the ions is used. The separation of this potential 
into a long-range electrostatic term and a short-range 
hydration term allows each of these two interactions to 
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be treated in a simple though approximated form. The 
virial expansion is a standard choice for approximating 
short-range interactions. Such an expansion fails for the 
electrostatic interaction due to its long-range . On the 
other hand, the wide success of PB theory demonstrates 
that the electrostatic interaction can be treated quite well 
in the mean-field approximation. Therefore we use this 
approximation for the long-range part of the interaction, 
and in this respect we remain within the framework of 
PB theory. 

The free energy (|l|) can also be obtained by rewriting 
the grand canonical partition function as a field-theory 
partition function. The short-range part of the interac- 
tion can be separated from the electrostatic interaction 
and a different expansion can be performed for each of 
these interactions. By using a density expansion for the 
short-range interaction and a loop expansion for the elec- 
trostatic interaction, Eq. (|^) is obtained up to second 
order in the density expansion and first order in the elec- 
trostatic potential [Q. 

The simplicity of our approach can lead to elegant an- 
alytical results, but has several limitations. The use of 
only the second term in the virial expansion implies that 
we are using a low density approximation. The valid- 
ity of such an approximation for a bulk fluid can be as- 
sessed by considering kBTB2C, where B2 is the second 
virial coefficient in the expansion of the pressure, and 
c is the ion density. Qualitatively, if B2C is small com- 
pared to unity, the correction to the ideal gas behavior is 
small and truncating the virial expansion after the second 
term is sensible. For non-homogeneous cases, the corre- 
sponding quantity is (1/2) J2j J dr' c(r')J7y(r — r'). For 
relatively high surface charges considered here this inte- 
gral approaches values of order unity near the charged 
surfaces, indicating that the approximation should only 
be expected to give qualitative results. Another defi- 
ciency of the virial expansion to second order can be 
seen from the fact that the direct correlation function 
is simply —Uij(r). This implies that the hard core inter- 
action is not described accurately in our treatment. A 
faithful description would require the vanishing of the 
pair correlation function /i2(r) for separations smaller 
than the hard-core diameter. Hence, it should be kept 
in mind that our main concern is to study the effects 
of a short-range interaction with a dominant attractive 
part. Finally, the fact that we describe the electrostatic 
interaction in the mean-field approximation implies that 
ion-ion correlations are ignored, as they are in PB theory. 
When our approach is applied for the aqueous pair po- 
tential model, these approximations should also be kept 
in mind. In particular, we follow Ref. |34| and do not 
include an effective ion-surface potential |[49|| . 



B. Density equations 

The mean density distribution is obtained by minimiz- 
ing the total free energy fl = flpB + Ail. From equations 
(I), (^ and (|l|) we have: 

^ = ^1 (V*)'d3r + kBTjj2 ^'^ (log f - l) d'r 
+ ^E/^^W^'.('--0^.('-')d^rdV 

+ j A(r) (^H + T ^ ""'^^ ^^^^ 

where and Q — exp(/3/ii)/Ai^ are the chemical po- 
tential and the fugacity of the ion species i, respec- 
tively. The thermal de Broglie wavelength, At, is equal 
to h/ , where h is the Planck constant and 
m is the ion mass. Requiring an extremum of il with 
respect to ^ gives: A = (e/47r)^' as in Eq. ([Tl|). Taking 
the variation with respect to Ci then gives: 

log ^ + E / c,(r')[/y (r - r') d\' + Pe,^{v) = 

(20) 

This equation is supplemented by the Poisson equa- 
tion (|^). Since Eq. ( po|) is an integral equation, the q 
cannot be written as a simple function of \E' as in the 
PB case. Therefore, a single equation for analogous 
to the PB equation, cannot be obtained, and we are left 
with the two coupled integro-differential equations ( |20| ) 
and (0). These equations should be solved together to 
obtain the electrostat ic p otential and density profiles. In 
the case J7 ^ 0, Eq. ( pof ) reduces to the Boltzmann rela- 
tion Ci = Ci exp(— /Jci^*) with Qi = Ch.i- Combining this 
relation with Eq. @) reproduces the PB equation (p^. 

In order to simplify the set of equations, we assume the 
same short-range interaction between the different pairs 
of ion species. Assuming that the charged surfaces are 
negatively charged, we choose: Uy(r) = u+_|_(r) = u{r), 
where U4._|_(r) is the short-range effective potential be- 
tween the (positive) counter-ions. This assumption is 
not exact for the effective potentials of ions in water . 
However, since only the counter-ions reach high densities, 
close to the oppositely charged surfaces, and the co-ions 
are repelled from the surface neighborhood, the exact 

choice of the potentials (r) and u (r) is expected 

to be of only minor significance. 

We now consider an electrolyte of valency z+:z_, i.e., 
a solution of positive and negative ions of charges e± = 
±z±e, where e is the electron charge. We designate 
the surface charge density on the plate as a constant 
a and the bulk densities of the positive and negative 
ions as Cb = Cb,+ and Cb,-, respectively. Due to charge 
neutrality in the bulk, Cb,_ = (z+/z_)cb and similarly, 



6 



C_ = (z_(-/z_)C where C = C+- Equation ( po| ) can then 
be written as foUows: 



c±(r) = C±e^''"'±*exp 



J c{r')U{r-r')<fr' 



(21) 



where c(r) = c+(r) + c_(r) is the total ion density, and 
U{r) = U++{r) is obtained from u(r) using Eq. (17). 
From the Poisson equation (^ we obtain: 



(Z+C4 



47re 



Cz+ (e 



X exp 



-J c(r')C/(r-r')d3r' 



(22) 



Note that in addition to the exphcit dependence on the 
ion valencies z± in equations (|l|) and (^), m a more re- 
alistic model the details of the potential u{r) should also 
depend on the type of counter-ion species present in the 
problem. 



III. SINGLE CHARGED PLATE 
A. Density equations 

After presenting the general formalism let us consider, 
as an example, a single negatively charged planar sur- 
face (Fig. 1). The charged surface is in contact with an 
electrolyte of valency z+:z_. We designate the axis per- 
pendicular to the plate as the z axis, and consider the ion 
solution in the region z > 0. For simplicity we consider 
positive and negative ions of the same hard-core diame- 
ter dhc- The coordinate of closest approach of the ions to 
the plate is designated as z = 0. Hence the "real" surface 
lies at a distance of one ion radius cihc/2 from the actual 
z = plate position, as shown in Fig. la. When we refer 
to conventional PB results, however, the ions are point- 
like and the plate should be understood to be positioned 
exactly at z = 0. 

Due to the one-dimensional symmetry imposed by the 
uniformly charged planar plate, the integration in Eq. 
( pT| ) can be performed over the x — y plane, leaving us 
with profiles depending only on z, the distance from the 
plate: 



c±{z) = C±e 



exp 



c{z')B{z - z')dz' 



(23) 



where c(z) = c+(z) + c_(z) is the total ion density and 
B{z) is the integral of U{r) in the plane of constant z. 
Using cylindrical coordinates: 



B{z) = 2tt f pdpU (^^/z^ + p^^ 



(24) 



and the Poisson equation (p3) reads: 



d^ 

dz2 



Ane 

e 



Cz+ (e'^'^^-* 



X exp 



c{z')B{z - z')dz' 



(25) 



Equations (|2^) and (^5|) are supplemented by the bound- 
ary conditions: 



d* 

dz 



47r 



z=0 



d7 







(26) 



Finally, the relation between C and the bulk density 
Cb can be obtained from Eq. (23). As z oo, be- 
comes zero, and c± assume their asymptotic constant, 
bulk values. Thus the integrand inside the exponential 
can be replaced by —(1 -I- z+/z_)cb-B(z — z'). Recalling 
that C-i- = Cb and c_ = (z-|_/z_)cb, we obtain: 



Cb = C exp 



z_ 



BtCh 



where: 



Bt= I dzB{z) = / d?YU{v) 

J —OO J 



(27) 



(28) 



is also equal to 2_B2, the second virial coefficient. Note 
that B(z) and i?t are negative for an attractive inter- 
action. The limit -BtCb — > is the limit in which the 
short-range interaction becomes negligible in the bulk. 
In this limit the relation between the bulk density and 
fugacity of Eq. ( p^ ) tends to the ideal gas relation Cb = 
C = exp(/3Af)/A|. 

Two special cases will be of particular interest in the 
following sections. The first is the case of a monovalent 
1:1 electrolyte, where we have: 



c±(z) = Ce^^^* 



exp 



c(z')S(z - z')dz' 



d^ 

dz^ 



47r 

-c(z) 



and: 



Cb = Cexp(-2BtCb) 



(29) 



(30) 



The second case is that of no added salt. The solution 
contains only monovalent counter- ions (z+ — 1, z_ = 0). 
This case can be obtained by taking formally the limit 
C — > of Eq. (^9|), or by repeating the derivation from 
Eq. ([l9| ) with only one type of ions, of charge e. The 
term —k^T J d'^rclog(C) in fl is then a Lagrange multi- 
plier added to impose the condition: dzec{z) = \a\. 
The following equations are then obtained: 



c(z) = Coe-'^'^* 



exp 



c{z')B{z — z') dz' 



d^* 47re , , 
di^ ^ -—'^'^ 



(31) 
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where Co is an arbitrary reference fugacity. The choice 
of Co determines the (arbitrary) position in which 5* is 
zero. Note that the electrostatic potential 4" diverges in 
the bulk. This divergence exists also in the usual PB 
theory, because the system is effectively one dimensional 
with no screening by added salt. Although ^'(z) has a 
weak logarithmic divergence, the density of counter-ions 
decays to zero, lim^^oo c{z) = as it should. 



B. Parameters and length scales 

For the ion- ion potential u{r ~ r') we use an effective 
potential between Na"*" - Na"*" ion pairs. The potential 
was calculated using a Monte-Carlo simulation ^ for 
an NaCl ionic solution of concentration 0.55 M, at room 
temperature. The electrostatic interaction between the 
ions is subtracted, and the net short-range potential is 
shown in Fig. 2. For ion-ion separations below 2.9 A a 
hard core interaction is assumed. Fig. 3 shows the func- 
tion B{z), derived from this potential, using Eq. (24). 
Note that B{z) has several local maxima and minima. 
These correspond to the local maxima and minima of 
u{r). Thus the structure of B{z) reflects the oscillatory 
behavior of the effective potential. 
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FIG. 3. The effective interaction in a planar geometry B{z) 
obtained from the potential of Fig. 2, using Eq. (^^. The os- 
cillating structure of the radial potential shown in Fig. 2 is 
apparent in the secondary minima of B{z). 

We use the effective potential calculated for Cb = 
0.55 M, regardless of the actual bulk ion concentration 
in the system. Since the important effects occur near 
the charged surface, where the ion concentration is much 
larger than Cb, it seems reasonable to use an effective 
potential calculated in the presence of a rather high salt 
concentration. The choice of Cb = 0.55 M is still some- 
what arbitrary, and we rely on the fact that the depen- 
dence of the effective potential on the ion concentration 
is weak 

It is useful to employ two length scales that charac- 



terize the PB density profiles g. The Gouy- Chapman 
length, defined as 6 = ek-BT/{2'!Te\a\), characterizes the 
width of the diffusive counter-ion layer close to a single 
plate charged with a surface charge cr, in the absence of 
added salt. The Debye-Hiickel screening length, Ad = 
(87^CbeV£fcB^)"^/^ equal to 19.6 A for Cb = 0.025M at 
room temperature characterizes the decay of the screened 
electrostatic interaction in a solution with added salt. 
The strength of the electrostatic interaction can also be 
expressed using the Bjerrum length. Is = e^Kek-QT). 
This is the distance at which the electrostatic interaction 
between two unit charges becomes equal to the thermal 
energy. The Bjerrum length is equal to about 7 A in wa- 
ter at room temperature. 

The inclusion of hydration interactions introduces ad- 
ditional length scales in the system. For the interaction 
shown in Figs. 2 and 3, the range of the interaction dhyd 
can be seen to be approximately 7 A, over twice the hard 
core diameter dhc = 2.9 A. The strength of the hydra- 
tion interaction is characterized by Bt ~ —(7.9 A)" 
calculated from Eq. 



C. Numerical results 

Equations (^3|) and (|2^) are a set of three nonlinear 
integro-differential equations. We treat them numeri- 
cally using an iterative scheme, based on the assump- 
tion that the positive ion density profile is dominated by 
the electrostatic interaction. We start with the analyti- 
cally known PB profile close to a single charged plate and 
calculate iteratively corrections to this profile, as result 
from equations ( p3| ) and (p5|). For a 1:1 electrolyte we 
iteratively solve the equation: 



d2vl/(") 



Bttc 

e 



Csinh {pe¥"^^ 

c^''-^\z')B{z-z')dz' 



X exp 



(32) 



where c(z) = c+(z) -|- C-{z) is the total ion density and 
the superscript n stands for the nth iteration. For n > 0: 



X exp 



=("-i)(z')S(z-z')dz' 



(33) 



and the zeroth order densities c^'' are taken as the den- 
sity profiles generated by the PB equation (|l^). The 
boundary conditions (|2^) are satisfied by the electrostatic 
potential 5''^"-' in all the iterations. Note that using our 
iterative scheme, Eq. (^2|) is an inhomogeneous differen- 
tial equation, because the integral in the exponential is a 
known function of z, calculated numerically in the (n — l) 
iteration. A similar iterative scheme, based on Eq. ( |3l| ) 
can be used when only counter-ions are present in the 
solution. 
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FIG. 4. Counter-ion density profile (solid line) obtained 
from numerical solution of Eq. ( |3l| ) with the hydration in- 
teraction as of Fig. 3, plotted on a semi-log plot. No 
salt is present in the solution. The surface charge is 
\a\ = 0.333 C/m^. The dielectric constant is e = 78 and 
the temperature is T = 298 K. The dotted line shows the 
corresponding density profile obtained from the PB equation. 

Figure 4 shows the calculated density profile of the 
counter-ions on a semi-logarithmic scale, for a charged 
plate with a surface charge, \a\ = 0.333 C/m^, corre- 

sponding to an area of approximately 48 A per unit 
charge. This is a typical high surface charge obtained 
with mica plates. It corresponds to a Gouy-Chapman 
length b — 1.06 A, at a temperature of 298K, with e = 78. 
No salt is present in the solution. The calculated den- 
sity profile (solid line) is compared to the PB prediction 
(dotted line). The short-range attraction favors an in- 
creased concentration of counter-ions in the vicinity of 
the charged plate. This results in an increase of the con- 
centration relative to the PB prediction. For a surface 
charge as in Fig. 4, an increase of the concentration is 
seen at distances from the plate up to approximately 
4.5 A. The overall number of counter-ions is fixed by the 
requirement of charge neutrality. Therefore, the increase 
in the density of counter-ions near the plate is balanced 
by a reduced concentration further away. 

When salt is present in the solution, the short-range 
attraction draws additional ions from the bulk solution 
to the diffuse electrical layer near the plate. This can 
be seen in Fig. 5, in a comparison of counter- ion pro- 
files for different values of the bulk concentration Cb- For 
each salt concentration, the figure shows the ratio be- 
tween the counter-ion density and the density predicted 
by PB theory, as a function of the distance from the plate. 
The dotted line shows the result in the no-salt limit. As 
the salt concentration increases, the counterion concen- 
tration increases relative to the PB concentration at all 
distances from the charged plate. Qualitatively, however, 
the hydration effect on the counter-ion profile is similar 
in all the curves. As long as the Debye-Hiickel screening 
length is large compared to the Gouy-Chapman length. 
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FIG. 5. The ratio c+/c+^ between the positive ion density 
obtained from Eq. ( p^ ) and the value obtained from PB the- 
ory, for a surface charge \a\ = 0.333 C/m^ and several values 
of Cb. Other parameters are as in Fig. 4. The three values of 
Cb: 0.1 M, 0.025 M and 0.00625 M correspond to Debye-Hiickel 
screening lengths Ad — 9.8 A, 19.6 A and 39.2 A, respectively. 

b — 1.06 A, the density profile in the vicinity of the plate 
is dominated by the balancing counter-ions and the salt 
has only a small effect. 
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FIG. 6. The ratio c+ /c^ between the positive ion den- 
sity obtained from Eq. (E3) and the value obtained from PB 
theory, for surface charges \a\ — 0.333 C/m^ (dashed line), 
0.1 C/m^ (solid line) and 0.0333 C/m^ (dotted line). The bulk 
salt concentration Cb is 0.025 M. Other parameters are as in 
Fig. 4. 

The effect of the hydration interaction is strongly de- 
pendent on the surface charge a. As a is increased, the 
ion density near the surface increases too. The exponen- 
tial in Eq. (^ij) deviates more strongly from unity, leading 
to a larger deviation from PB theory. The dependence 
on cr is demonstrated in Fig. 6. The ratio of the positive 
ion density to its PB value is shown for three values of 
the surface charge. The effect of the hydration potential 
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is very minor for small surface charge, |cr| = 0.0333 C/m^ 
(dotted line), where the deviation from PB is less than 2% 
at its maximum, and considerable for a surface charge of 
0.333 C/m^ (dashed line), where the deviation from PB 
reaches almost 40%. 

The numerical scheme, described above, requires sev- 
eral iterations to converge fully. It is interesting to note, 
though, that the first iteration captures most of the ef- 
fect of the short range interaction. This indicates that the 
density profile is dominated, as we assumed, by the elec- 
trostatic interaction, and assures that the convergence of 
the iterative scheme is good with the PB density profile 
as the zero-th order approximation. On the theoretical 
level it indicates that the effect of the hydration interac- 
tion can be seen as a perturbation over the PB results. 
The fact that the first iteration provides a good approx- 
imation to the full iterative result can lead to further 
analytical approximations. For example, the corrections 
to the density profiles, in the no added salt limit, are 
studied analytically in the next section, based on this 
observation. 

As an example for the results of the first iteration, we 
compare, in Fig. 7, the correction to the counter-ion den- 
sity profile obtained in the first iteration (dashed line), 
with the full iterative result (solid line). We use a high 
surface charge of 0.333 C/m^, where the differences be- 
tween the exact profile and that of the first iteration are 
relatively pronounced. The two density profiles differ by 
at most 3.2 percent, where the ion density deviates from 
the PB value by 30 percent. For smaller surface charge 
the results obtained in the first iteration are even better. 




FIG. 7. The positive ion density profile obtained after one 
iteration of E g. ( ^2[) (dotted line), compared to the full so- 
lution of Eq. (E9|) (solid line). Parameters are as in Fig. 5. 
The maximal deviation between the two density profiles is 
3.2 percent, where the deviation from PB is approximately 
30 percent. 



D. Contact density and the contact theorem 

The contact density of the ions is barely modified as 
compared with the PB prediction. This is evident in 
Figs. 4-6. As long as the Debye-Hiickel screening length 
is large compared to the Gouy-Chapman length, or the 
hydration interaction is negligible in the bulk, the modi- 
fication remains small. This result can be obtained from 
a generalization of the PB contact theorem ||,^^ : 



27r/3 2 

e 



PhuU 



(34) 



where Pbuik is the bulk pressure of the ionic solution. 
Equation (^^ is derived in detail for the free energy 
used in our model in Ref. It is obtained from the 

equality of the internal pressure in the electrolyte solution 
at different distances from the charged plate. Far away 
from the charged plate the pressure must be equal to the 
bulk pressure of the ionic solution, because the densities 
approach their bulk values and the electrostatic poten- 
tial becomes constant. At the contact plane between the 
plate and the solution, the pressure involves only an elec- 
trostatic contribution and an osmotic contribution, as in 
PB theory. This is due to the fact that in our model no 
short range interaction between the plate and the ions is 
included. Equating the pressure at the contact plane and 
far away from the plate results in Eq. ( p^ ) . 

The contact density, as expressed by Eq. (|3^), differs 
from the PB prediction only due to the change in the ac- 
tual value of Pbuik- This change is negligible if the short- 
range interaction is not of importance in the bulk. In 
addition, if the surface charge is high, such that 6 <C Ad , 
Pbuik is negligible compared to the second term in the left 
hand side of Eq. (|3j) . Thus the contact density remains 
very close to the PB prediction. In the no-salt limit Pbuik 
is zero and the contact density coincides exactly with the 
PB result, c+(0) = {2TTf3/e)(T^. 



IV. ANALYTICAL SOLUTIONS 

The simplicity of the model makes it possible to obtain 
various analytical results. The effect of the hydration on 
the ion distribution can be characterized by several quan- 
tities, such as the magnitude of the deviation from the PB 
result and the effective PB surface charge density seen at 
a distance from the plate. Using several simplifying as- 
sumptions it is possible to obtain analytical expressions 
for these quantities. 

First we assume that the hydration interactions can 
be neglected in the bulk, i.e., PtCb <C 1. In this case, 
the effect of the hydration potential is significant only in 
the vicinity of the charged surface, where the ion density 
becomes large. In addition, the Debye-Hiickel screening 
length. Ad, is taken to be large compared to the Gouy- 
Chapman length b = e/(27rZB|cr|). Since Ad ^ b, the 
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negative co-ion density near the negatively charged sur- 
face can be neglected compared to the positive counter- 
ion density. Far away from the charged plate, the system 
is well described using the PB equation, with an effec- 
tive surface charge density dcff different from the actual 
charge density cr. The result of the above two simplifying 
assumptions is that the salt is of minor importance in the 
region where the effective surface charge is determined. 
The effective surface charge can then be inferred by con- 
sidering the case in which only counter-ions are present 
in the solution (no added salt). 

Equation (|l]) can now be recast in a simpler form, by 
considering r] = log(c/Co), as expressed by Eq. and 
taking its second derivative: 



^° dz2 



The PB density profile, cpb(z) = Cqc''"^''^, for the 
same surface charge, satisfies the equation cPrjo/dz'^ = 
{4:7rPe'^(o/e) exp{j]o). Its exact solution is known to be: 



cpBiz) = Coc 



1 



27r/B {z + by 



(36) 



Note that only in the PB equation ri{z) is the reduced 
electrostatic potential e^'(z)/fcBT. From the general- 
ized contact theorem (p^), the surface density in the 
no added salt case and in the presence of one plate is 
c(0) = 2TrPa^/e, as in PB theory. Therefore: 



■i]{z = 0) = rio{z = 0) 
From the derivative of c(z), Eq. (^l|), we find: 



(37) 



This equation can be further simplified by omitting w{z') 
from the integrand in the right hand side. This approx- 
imation was motivated in Sec. |III C| and is equivalent to 
stopping the iterative scheme (|32|) after the first itera- 
tion. The density profile is then replaced by the PB 
density profile in the term that involves the hydration 
interaction B{z). This results in the equation: 



(ie^cpB{z)w{z)+V{z) = 



(42) 



where T{z) is the convolution integral 
1 



r(z) 



dz' (43) 



27r/B Jo {z' + b) 



dz2 



The corresponding boundary conditions, obtained from 

same approximations. 



I ine corresponding boundary com 

dz (35) equations and (|3|) using the 



are: 



w{z 
dw 
d7 



0) 



, „ dB{z') 

dz cpb(z ) — -, 

dz 



(44) 



Equation (|4j) is a second order linear differential equa- 
tion for w{z) and can be solved analytically. The solu- 
tion, given in detail in Appendix is expressed in terms 
of the convolution integral r(z) of Eq. (43). The effec- 
tive surface charge and the effect of the hydration on the 
density profile can then be calculated in several limits, 
described in detail in Appendix Here we outline the 
main results. 



A. Slowly varying density: b 3> dhyd 



dry d* 
dz dz 



, ,, dB{z-z') 

dz dz ) 

dz 



and using the boundary condition (Eq): 



dry 




dryo 




J 


dz 


2 = 


dz 


4 

z=0 


Jo 



' , dB(zO 

dz c(z ) — 

dz 



(38) 



(39) 



where the odd parity of d_B/dz has been used. This re- 
lation can be used together with Eq. ( |37| ) as a second 
boundary condition at z = 0, instead of the boundary 
condition of vanishing dry/dz at infinity. 
Linearizing Eq. (p5[) with respect to: 



(40) 



w = 7] -ria = log(c/cpB) 



which is valid for relatively small deviations from the PB 
profile, results in the following equation: 



d w 47r 2 / N / N 

-TT CPB[Z)W[Z) 

dz'' e 



d^Blz - z') 
dz'{l + w{z'))cMz') ' 



(41) 



In the limit h » rfhyd, the PB distribution varies 
slowly on the scale of the hydration interaction, de- 
scribed by B{z), and the theory becomes effectively a 
local density functional theory. The specific form of 
B{z) is not important, and all the results simply de- 
pend on Bt = B{z) dz. The deviation of the ef- 
fective Gouy-Chapman length bcs from the actual Gouy- 
Ghapman length b depends linearly on Bt and on the 
surface charge a ^ 1/b. This can be expected since we 
use a linearized equation. Thus we have, on dimensional 
grounds, bcs — b ^ Bt/lBb- The detailed calculation gives 
the numerical prefactor: 



-Bt 1 



Oeff 



47rZB b 



(45) 



Since Bt is negative bcs is larger than b and the effective 
surface charge, Uoff, is smaller than the actual surface 
charge a. This result should be expected. The short 
range interaction attracts counterions to the vicinity of 
the charged plate and the surface charge is screened more 
effectively than in the PB equation. 

The correction to the counter-ion density profile, de- 
scribed by w{z) = log[c(z)/cpB(2;)], is found to be: 
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w{z) 



1 



-Bt 



2ttIb {2{z + by b{z + b) 



(46) 



The density profile is increased relative to PB theory 
for distances smaller than b/2, and decreased for larger 
distances. The deviation from PB, wlz), is maximal at 
z = 0, where it is equal to —Bt/^AnlBb^), and minimal 
at z = 2b, where it is equal to Bt / {12nl'Qb^) ■ 

Figure 8 shows the approximated function 'w{z) of Eq. 
( p6| ) for b = 21.2 A, corresponding to &/dhyd ~ 3 (dotted 
line). The approximation is compared with the function 
w{z) obtained from the exact solution of equation ( |3l| ) 
for the case of no added salt (solid line). Although b is 
not much larger than dhyd, the approximation describes 
well the correction to the PB profile. Note that 'w{z), 
as expressed by Eq. (|4^) is maximal at z = 0, whereas 
according to the contact theorem w{Q) should be zero. 
This apparent inconsistency results from neglecting the 
range of the hydration potential relative to b. In the pre- 
cise solution of E q. (pi^ ) ?«(0) is zero, as it should be. The 



prediction of Eq. (|46|) is valid only for distances z > dhyd, 
as can be seen in Figure 8. 

The range of validity of the linearization procedure can 
be found by requiring that the minimal and maximal val- 
ues of w{z) are small compared to unity: 



-St 



< 1 



(47) 



scale of the hydration interaction. The effective Gouy- 
Chapman length has the same form as in the limit of 
slowly varying density, b 3> rfhyd, but having a different 
prefactor: 



Ooff 



-Bt 1 
127r/R b 



(48) 



The effective surface charge is, therefore, smaller than the 
actual surface charge. Note that bcs depends on B{z), 
in this limit, only through B^. The linear dependence on 
cr ^ 1/5 follows from the linearization leading to Eq. (|4^), 
as in the previous limit. 

It should be stressed that although b is small compared 
to dhyd we still assume that b is large enough for the lin- 
earization to be valid, i.e., we assume that w{z) is small 
compared to unity. Furthermore, the counter-ion density 
should be small enough that we can sensibly use only the 
quadratic term in the virial expansion. To check the va- 
lidity of these assumptions, the correction to the density 
profile should be considered. 
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FIG. 8. The logarithm of the ratio between the counter-ion 
density obtained with the inclusion of the hydration interac- 
tion and its value in PB theory, w{z), as a function of the 
distance from a charged plate, with no added salt in the so- 
lution. The solid line shows the function w{z) obtained from 
the exact solution, for b = 21.2 A. The dotted line shows the 
approximated curve obtained from the linearization with re- 
spect to Eq. (Ea), in the limit h S> dhyd, Eq. (Eg). 



B. Surface layer limit: h <C dhyd 

In the limit in which b <^ rfhyd, the ion density effec- 
tively becomes a dense layer concentrated at z = on the 



The form of w{z) depends, in the surface layer limit, on 
the specific form of B{z). In order to study w{z) analyti- 
cally, we use an approximated form of B{z), described in 
Appendix A typical form of the approximated w{z), 
obtained using this approximation [Eq. (B12)], is shown 
in Fig. 9 (dotted line). The Gouy-Chapman length is 
b = 1.06 A, corresponding to 6/(ihyd ~ 0.15. In addition, 
the function w{z) obtained from the exact solution of 
equation ( pl| ) is shown for comparison (solid line). The 
approximated curve captures well the qualitative behav- 
ior of the correction to the PB profile. Note that the 
discrepancy between the approximated and actual pro- 
files results not only from the linearization and small b 
limit, but also from the loss of detail due to the use of an 
approximated form for B{z). 

The deviation from the PB profile, w{z), can be quali- 
tatively described as follows. For z < dhc, w{z) increases 
from zero quadratically (with an additional term of the 
form z^logz) to its value at z = rfhc- It then decreases 
from its maximum positive value to a minimum, nega- 
tive value, on a scale of the range of the attractive part 
of B{z). This minimum value is equal to approximately 
Bt /6'KlBbdhyd- For distances larger than the interaction 
range, w{z) assumes the form w{z) ~ 1/z, characteriz- 
ing a PB profile with a modified, effective surface charge. 
For finite values of 6, we can expect the above behavior 
to be smoothed over a scale of order b. 
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FIG. 9. The logarithm of the ratio between the counter-ion 
density obtained with the inclusion of the hydration interac- 
tion and its value in PB theory, w{z), as a function of the 
distance from a charged plate, with no added salt in the so- 
lution. The solid line shows the function w{z) obtained from 
the exact solution, for b = 1.06 A. The dotted line shows the 
approximated curve obtained from the lineariza tion with re- 
spect to w, Eq. (^^, in the limit b <C dhyd, Eq. ( B12 ). 



The validity of the linearization can be found by re- 
quiring that 1^(2)1 <C 1. This requirement resuhs in the 
following condition: 



< 1 



(49) 



The vahdity of stopping the virial expansion at the 
quadratic order can be shown to have the same condi- 
tion. For the hydration potential of Fig. 2, the condition 
expressed in Eq. ( p9| ) implies that the various approxima- 
tions we use start to break down when b becomes smaller 
than approximately 1 A, or a > 0.022 e/A . When & is of 



this order, it is well below dhyd, making the surface layer 
limit a sensible approximation. 



C. Effective surface charge 

In the two limits described above, the effective Gouy- 
Chapman length was found to be of the form 6off — b ^ 
—Bf /lBb, with different prefactors in the two limits. For 
intermediate values of b, the effective charge depends on 
the specific structure of the function B{z). In order to 
study this dependence, we use a simple approximated 
form for B{z), described in Appendix ^ Using this ap- 
proximation, an analytical expression can be obtained for 
the effective Gouy-Chapman length for all values of b. 

Figures 10a and 10b show the predicted bcs and bcS — b, 
respectively (solid lines) as a function of &, together with 
the asymptotic limits ( ^ ) and (48) (dotted lines). As 
the surface charge increases from zero (and b decreases 
from infinity), the effective charge |croff| increases too 
(but is always smaller than the actual surface charge). 
When b reaches a certain value 6""", bcs starts increas- 
ing with further reduction of 6, i.e., the effective charge 
decreases with increasing surface charge above |(t|™^'^ = 



e/(27rZB6™'"). The value of 6™™ depends on the structure 
of the function B{z), but can be estimated to be between 
the values predicted by the asymptotic expressions ( |45| ) 
and (Eq). From the condition d&cfi/d6|b^bmin = we find: 



and: 



-B, 



< b" 



< 



-Bt 
47r?R 



"eff 



2b" 



(50) 



(51) 
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FIG. 10. The effective Gouy-Chapman length bctt (a) and Ab = b^s — b (b), as a function of the Gouy-Chapman length b. 
The solid lines show the behavior predicted by Eq. (p314[ ), with Bt = —500 A , dhc = 2.9 A and Bo = 41.8 A . The dotted lines 
show the asymptotic limits of equations (^) and (|^. The symbols show results extracted from numerical solutions of Eq. 
(0), using B{z) of Fig. 3, with salt concentrations of 10~^ M (circles) and 0.1 M (crosses). The salt has a very small effect. 
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For the hydration interaction of Fig. 2, Bt is approx- 
imately —500 A . The value of 6™™ is then between 
1.36 A and 2.35 A, corresponding to a surface charge den- 
sity between 0.15 C/m^ and 0.26 C/m^. The values ob- 
tained from the approximated curve, shown in Fig. 10, 
are 6™ ~ 1.5 A and 6™" ~ 3.4 A. 

For small enough values of b, the effective surface 
charge |crefi| should increase again with an increase of 
|(t| and become larger than \a\. This effect cannot be 
predicted by our model because of the low density ap- 
proximation used for the hard core interaction. In par- 
ticular, the hard core of the ions should cause the den- 
sity to saturate at the close packing density, leading to 
a reduced screening of the surface charge relative to PB 
theory |2^j2^,|5l|. In our model, as in the PB theory, 
the counterion density near the surface is not bounded, 
and increases indefinitely as a is increased. Although our 
model includes the steric repulsion between ions, this re- 
pulsion is "softened" , and is always outweighed by the 
attractive part of the ion-ion interaction. 

In addition to the prediction obtained using the lin- 
earized approximation. Figure 10 shows values of bcs 
extracted from numerical solutions of the full equation 
(p9|), using the original interaction B{z). The equation 
was solved with two different salt concentrations: 10^^ M 
(circles) and 0.1 M (crosses). The value of 6eff was esti- 
mated from the positive ion density at large distances 
from the plate, by finding the value of h that would re- 
sult in the same calculated values of the density in a 
solution of the PB equation. Note that for both salt 
concentrations, 6cff is very close to its predicted value, 
meaning that the salt has a very small effect on Ucff. 
This result is not obvious for the high salt concentration 
of 0.1 M. The Debye-Hiickel screening length is approxi- 
mately 9.6 A, not much larger than the range of the hy- 
dration interaction, dhyd — 7 A, and comparable to the 
Gouy-Chapman length at the large h region of the plot. 



V. CONCLUSIONS AND OUTLOOK 

In this work we have studied the effects due to the dis- 
creteness of the solvent in aqueous ionic solutions. Hy- 
dration interactions are found to have a significant ef- 
fect on the structure of the diffusive layer near highly 
charged surfaces. The counter-ion density is increased in 
the vicinity of the charged surface, relative to the PB pre- 
diction, and decreased further away. The distance from 
the charged plate in which the density is increased, and 
the magnitude of the deviation from the PB density, de- 
pend strongly on the surface charge, and on the parame- 
ters of the short-range hydration interaction between ion 
pairs. 

The ion-ion hydration interaction can be described 
roughly using two parameters. The first parameter is the 
range of the hydration interaction, c?hyd, equal to approxi- 
mately 7 A for Na+-Na+ pairs. The second parameter, Bt 



has dimensions of volume and characterizes the strength 
of the hydration interaction. It is equal to approximately 
-500 A for Na+-Na+ pairs. Two limits can be consid- 
ered, where the Gouy-Chapman length, h ~ l/cr, is small 
or large compared to the range of the hydration inter- 
action dhyd- In both of these limits we assume that the 
Debye-Hiickel screening length. Ad , is large compared to 
h and dhyd- 

In the limit b 3> rfhyd, the counter-ion density be- 
comes depleted, relative to the PB prediction, starting 
at a distance z ~ h/2 from the charged plate. The maxi- 
mum absolute value of w(z) — log[c(z)/cpB(z)] scales as 
—Bt/l-Qb^- In the limit h <^ dhyd, the distance from the 
plate, where the counter-ion density becomes lower than 
the PB prediction, is between z = c?hc and z = c?hyd- The 
maximum absolute value oi w{z) scales as —Bt/l^dhydb. 

Far away from the charged plate, the density profile 
can be well described using the PB theory with an ef- 
fective surface charge that can be calculated analytically. 
The correction to the Gouy-Chapman length in the two 
limits b ^ dhyd and b <C dhyd is always positive and 
scales as —Bt/l^b^ but has different numerical prefac- 
tors. When the surface charge on the plate is increased, 
the effective surface charge, dcff , is found to reach a cer- 
tain maximal value. Above this maximal value Ucs de- 
creases with further increase of the actual a on the plate. 
The various approximations we use start to break down 
when 6 is smaller than approximately — -Bt/67r^Brfhyd, cor- 
responding to 6 < 1 A. 

An important outcome of this work is that the correc- 
tion of the PB ion density due to the hydration inter- 
action is significant near highly charged surfaces. The 
electrostatic interaction dominates the ionic distribution 
and the hydration interaction can be seen as a perturba- 
tion. For a high surface charge density of, say, one unit 
charge per 48 A the counterion density deviates from its 
Poisson Boltzmann value by at most 30 percent. The 
effective change in the surface charge is more significant, 
from le/48A^ to about le/13A^. 

The hydration effect on inter-surface forces can be very 
pronounced, as opposed to the effect on the ion distribu- 
tion. This result will be presented elsewhere [Q. Our 
model predicts an attractive contribution to the pressure 
between two parallel charged plates. At distances below 
several nanometers this contribution can outweigh the 
electrostatic repulsion and lead to an overall attraction 
between the plates. Our two-plate findings can also be 
compared with available AHNC results showing 
good qualitative agreement both for the ion density pro- 
file and pressure. 

The formalism we present can be readily generalized 
to other geometries. This could lead to an estimation 
of the aqueous solvent effects on phenomena such as the 
Manning condensation on cylindrical polyions [52| , and 
charge renormalization of spherical mycelles or colloids 
1^. In this respect our formalism offers an advantage 
over the AHNC approximation which was applied so far 
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only in a planar geometry. Another interesting exten- 
sion of this work would be to consider the combination 
of fluctuation and hydration effects. This is particularly 
important for ionic solutions with divalent counter-ions, 
where fluctuation effects become large pl| , p^ , p^ . 
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APPENDIX A: INHOMOGENEOUS VIRIAL EXPANSION 



We consider an inhomogeneous system of particles with a short-range two-body interaction, and aim to express the 
free energy of the system in the low density limit as a functional of the density distribution. For simplicity we consider 
only one species of particles. The inhomogeneity of the system arises from the inclusion of an external field (fir), or 
from the boundary conditions imposed on the system. We begin by considering the grand canonical ensemble. The 
grand canonical partition function is: 




where fi is the chemical potential, At is the de Broglie thermal wavelength and Qn is: 

N 

Qjv= /nd'r,:e-'3^"«-^'» (A2) 



UnHv.}) = ifiin) + 1 ^ ^ « (|r, - r,|) (A3) 

i i j=ii 

We proceed on similar lines as the usual virial expansion in a bulk fluid, expanding log Zq in powers of the activity. 
Up to second order we have: 



+i(iJl) Jdhjd'r'c-l'('''>*^-'''»(e-''il'-''n-l) (A4) 

This can be seen as an expansion in powers of the fleld exp [/3 (/i — f{y^))] /A^- The local density c(r) can be expressed 
in a similar expansion: 

' " ' (A5) 



1 SlogZc 














-/3v(r) J 







This relation can be inverted to obtain an expansion of exp [/3 (/i — fir))] /X^ in powers of c(r). Up to the second 
order: 



g/3(M-v(r)) 



c(r) + c(r) f d^r' c(r') (l - e-^^d"^-""'!)) (A6) 



15 



and by substituting this relation in Eq. (A4) log Zq can be expressed as an expansion in c. Up to the second order: 

logZc = J d^rc{r) + ^ J d\ J d^r' c(r)c(r') (l - e^'^^d"^-'"'!)) (A7) 

The grand canonical potential can be obtained from the relation fl = — /ceTlogZc, with log^c given by Eq. (|A7|). 
In this expression, c(r) is the mean density profile for the imposed external field ip(j) and a given chemical potential 
fi. We would like to express SI as a functional of a general ion density c(r), whose minimization with respect to c(r) 
would give the equilibrium mean density. Regarding —k^T log Zq as a functional of x(r) = ip{r) — fj,, we have: 



SlogZc 
Sx{r) 



c(r) 



(A8) 



(A9) 



and expressing log Zq and x as functionals of c(r). We have already expressed logZc as a functional of c(r) in 
Eq. (A7). An expression for x(r) as a functional of c(r) can be obtained from Eq. (AE). Up to first order in c we 
have : 



The Legendre transform of this relation can be obtained by defining: 



e = -kBT log Zq - / d3rc(r)x(r) 



/3[(/.(r)-/i] = -log|A|c(r) 
= -log [A^cW] 



1 + J dVc(r')(l-e-Mk-'|)) I 
Vc(r') (l-e-Mk-r'Dj +0(c2) 



(AlO) 



Using this relation and Eq. (A7) we obtain, up to second order in c: 



/3e({c(r)}) = J d3rc(r) [log(A3c(r)) - l] 

r'c(r)c(r') ( 



+ i / d^r / d^r 



The functional O of c(r) has the property that 

se 



1 - e 



or equivalently: 



(5c(r) 



e 



d'^r c(r) [(p{r) — fi] 



m{c{r)}) 
Sc{r) 



(All) 



(A12) 



(A13) 



Thus, using Eq. (All), we obtain: 

r!({c(r)}) = fceT J dhc{r) (log^ ^ ^) + / 



-kjiT / d-^r / d-^r' 



r'c(r)c(r') ( 



■^r c(r)(^(r) 
1 _ e-^"(l''"'"'l) 



(A14) 



where C = exp(/3^)/A^. 

The derivation of Eq. ( A14 ) can be readily generalized to the case of several ion species of different charges and 
different pair interactions Uij (r) , resulting in Eq. ( p^ ) . 

A similar, more elaborate diagrammatic expansion of the thermodynamic potentials in the presence of an external 
field is presented in Ref. [Q. A variational principal for the grand canonical potential O is obtained in which Q is 
expressed as a functional of the mean density c(r) and the pair correlation function /i2(ri,r2). This expression is 
equivalent to Eq. (A14) up to the second order in the cluster expansion. 
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APPENDIX B: DETAILS OF ANALYTICAL RESULTS 



In this appendix we present details of the analytical approximations of Sec. IV. 

We consider first the analytical solution of Equation ( [42| ) . This equation is a second order linear differential equation 
for u'(z). Note that the function cpBiz) is a known function of z, given by Eq. (|36|). The solution of Eq. (42), with 
the boundary conditions of Eq. (H) is: 



w{z) 



1 



z + b Jo 



dz2 {z2 + b) 



dzi 



[zi + b) 



(Bl) 



where T[z) is the convolution integral, defined by Eq. By writing T{z) as: 



T{z) 



dz'T{z')5{z- z') 



(B2) 



w{z) can be rewritten in the following form: 

1 



w{z) = — 



z + b \ 3 
(z + b) 



b^r,,viz') 1 



6 3 



dz' {z' + bfT{z') 



z' + b 



(B3) 



The effective charge CTcff (or equivalently, the effective Gouy-Chapman length b^s) can be calculated from the coefficient 
of z^^ in w{z), as z approaches infinity: 



w{z) 



2 (b - 6eff) 



(B4) 



We thus find: 



bcs~b = - 

6 Jo 



dz 



&3 

Z + b 



-{z + bf 



r(z) 



(B5) 



A simple form for the convolution integral r(z) can be obtained in the limits in which b is small or large relative to 
dhyd, the characteristic range of the hydration potential. 



1. Slowly varying density: h ^ dhyd 



In the limit b ^ dhvd, the PB distribution varies slowly on the scale of the hydration interaction. The convolution 
integral V{z) of Eq. (g) can then be approximated in the following way: 



r(z) 



1 



27r/B 
1 

2^ 

Bt 



^ , H{z') d^B, 

dz TTTT . ^ \Z — Z ) 



277/ 



B 



dz' 

■OQ 

1 d6{z) 
62 dz 



(z' + 6)2 dz2 

1 d5{z 



62 dz 



6H{z) 



6H{z) 



(z + 6)4 



B{z 



(z + by 



(B6) 



wher e H {z) is the Heaviside function {H{z) — for z < and II(z) = 1 for z > 0). Inserti ng t his expr essi on in 
Eq. (B5) we obtain Eq. (^ ) for the effective Gouy-Chapman length. By substituting equation (B£) in Eq. (B3), the 
form of w(z), given in Eq. (fig), is obtained. 
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2. Approximated form for B{z) 



Some of the following results depend on the specific structure of the hydration interaction, characterized by the 
function B{z). In order to obtain analytical expressions, we use a simple approximated form, i3'*PP(z), instead of 
B{z). Assuming that the hydration interaction consists of a hard core interaction and a short-range attractive part, 
the function B{z) has some general characteristics that should be present in i3'^PP(z). For z < dhc, B{z) always has 
the parabolic form —{Bq + ttz^), where Bq — —B{z — 0). We assume that the attractive part of the interaction 
dominates over the short-range repulsion so that Bq is positive. For z larger than some finite value dhc + A, B{z) is 
practically zero due to the short range of the interaction. For < z < dhc + A, B{z) varies from —{Bq + Trd'f^^) to 
zero in a functional form that depends on the details of the attractive potential. The most simple way to model this 
behavior of B{z) is to have a linear increase of B^pp(z) between z = dhc and z = d^c + A, and to set B'^^p to be zero 
for z larger than dhc + A; 

r-(Bo + 7rz2) \z\<dhc 

B'-^nz) = - (So + TTdl) ^^^^+^-'^^ dhc < \z\ < dhc + A (B7) 
[o dhc + A<|2| 

The parameters in this expression should be chosen to match, approximately, the form of B{z). Setting dhc to be 
the hard core diameter of the real potential and setting Bq = —B{0) ensures that B(z) and B'^pp{z) are identical for 
z < dhc- The width A can then be set such that B^^^ ~ B^. 

2BQd^, + ^< + A (Bo + TTdl^) = -Bt (B8) 

This is desirable in light of equations and (Q), since the effective surface charge depends only on Bt in these 
limits. Figure 11 shows B{z) and i3''PP(z) for the hydration potential of Fig. 2. 



c<r- -20 




-80 h 

\ 



z[A] 

FIG. 11. The effective interaction in a planar geo metry, B{z), obtained from the potential of Fig. 2, and the corresponding 
approximated function £^^^(2), defined by Eq. (B7) (dashed line). The parabolic dependance for 1^1 < dhc is identical in the 
two curves. 



3. Surface layer limit: h <^ d 



hyd 



In the limit h <C dhyd, the convolution integral in Eq. ( |4^ ) becomes: 

|ct| d^B{z) 1 d?B{z) 



T{z) 



2'Klv,b dz 



The prefactor of T{z) in Eq. (B5) is + 0{b) and therefore the effective Gouy-Chapman length is: 



6 

feoff - b 



r dzz^B"{z)^^- 



(B9) 



(BIO) 
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This result is independent on the specific form of B{z). 

To obtain w{z), the deviation of the density profile relative to PB theory, Eq. (B9) can be substituted in Eq. 
Up to leading order in b the following expression is obtained: 



w{z) 



1 1 



dz' B"{z')z'' 



1 



dz'^B"{z') 



Using i3^PP(z), the approximated form of B{z) presented in the previous subsection, this gives: 



w{z) 



1 



GirlBb 
1 



47r 

Y 



"he 



3 ^ 



Bp + ndl 

rfhc(c^hc ^ 

Bo 



he 



- 27rlog 
A) z 



\z\ < dhc 



A(dhe + A)' 



dhc < \z\ < dhc + A 



Bt 1 



Qirl^b z ' 

The minimal, negative value of w{z) is assumed at z = dhc 

Bt 



J-hc 



A<\z\ 



w{dhc + A) 



A and is equal to: 
Bt 



67r/Bfc(dhc + A) 6nlBbd 



hyd 



(Bll) 



(B12) 



(B13) 



This results in the condition ( [49|) for the validity of the linearization in the surface layer limit. 

Using only the quadratic term in the virial expansion is sensible if /g dz' c{z')B{z — z') is small compared to unity. 
In the surface layer limit, this integral is simply: {\a\/e)B(z) = B{z) / {2'KlBb). Estimating the maximum value of 
(z)\ to be approximately — -Bt/(2dhyd) we obtain the requirement: — _Bt/(47rZB6(ihyd) ^ 1 , which is analogous to 

(I 



4. Effective Gouy-Chapman length 



Using B'^^^{z) in equations ( [43| ) and (B5) we find the following approximation for the effective Gouy-Chapman 
length: 



b.s-h= l-Bt^^ - ndlb + 27r4c62 + 2i?o log ( ^±^ili±^^) b 

b + dhc \ ^3 



(B14) 



— 27rlog 

This expression is shown in Fig. 10 and discussed in section p^V| . In the limits b iS> dhyd and b <^ dhyd it reduces to the 
asymptotic expressions (EH) and (48), respectively. 
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